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Abstract. Presented is a description of a Bayesian analysis framework for use with 
interferometric gravitational radiation data in searches for binary neutron star inspiral 
signals. Five parameters are investigated, and the information extracted from the 
data is illustrated and quantified. The posterior integration is carried out using 
Markov chain Monte Carlo (MCMC) methods. Implementation details include the use 
of importance resampling for improved convergence and informative priors reflecting 
the conditions expected for realistic measurements. An example is presented from 
an application using realistic, albeit fictitious, data. We expect that these parameter 
estimation techniques will prove useful at the end of a binary inspiral detection pipeline 
for interferometric detectors like LIGO or Virgo. 
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1. Introduction 



Great effort is currently being expended in the search for gravitational radiation. 
Detectors around the world pQ El El E] are reaching sensitivities where an event may 
be seen in the near future. The field of research has grown immensely since the first 
prediction of gravitational radiation j^j, or even from the time of the observations that 
confirmed its existence The inspiral of compact binary objects may provide the 
cleanest system for comparing observations with general relativistic predictions [7j. 
Observation of inspiral events could provide important information on the structure 
of neutron stars [HJ E]. Even cosmological information can be extracted from the 
observation of inspiral events fUl El El E3- The characteristics of radiation in the 
post-Newtonian regime will provide insight into highly non-linear general relativistic 
effects, such as the observation of the formation of a Kerr black hole as the binary 
system decays ^21E1E!- LIGO has recently conducted searches for signals from the 
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inspiral of binary neutron star systems in science run data [Tol 

Detection pipelines for binary inspiral signals CHI EH] m interferometer data 
will likely provide the first indication of detection. These pipelines provide parameter 
estimation through various estimation techniques. Other techniques, like those that we 
will describe in this paper, provide a means to examine the data and produce not just 
point estimates of parameters, but provide more detailed information e.g. on accuracy 
and correlations of parameters given the data at hand. The unique analysis methods 
that we bring to the examination of LIGO data are computationally intensive MCMC 
algorithms. Our method of parameter estimation is an exercise in Bayesian inference. 
These Bayesian statistical techniques offer great promise in problems where the number 
of parameters is large. Since its initial application in digital signal analysis 20 a MCMC 
methods have revolutionised many areas of applied statistics f2Tj . A distinct advantage 
of the MCMC approach is that computational time does not grow exponentially with 
parameter number, as it does for other methods [22]- After their introduction to 
the cosmological parameter estimation problem [2*3[ I24j . MCMC methods were quickly 
adopted and widely used with CMB data [23 1211 123, including WMAP [21 [23- We 
anticipate a similar scenario with gravitational radiation; after initial detection attention 
will be focused on parameter estimation and associated astrophysical information 
extracted from the data. The techniques that we describe in this paper are intended 
to be applied to interferometer data after an initial event has been registered by an 
inspiral detection pipeline our MCMC method would allow the extraction 

of information about the physical parameters. The technique we present here provides 
parameter estimates from the data output of a single interferometer. An MCMC 
technique using the coherent addition of signals from multiple detectors will be presented 
in a future publication. 

In the present paper we describe the MCMC routine we have developed to examine 
interferometer data for binary neutron star inspiral signals. We have implemented 
advanced MCMC methods in order to effectively and efficiently search a five-parameter 
space, find the signal's parameter values, and quantify and characterise them in detail. 
There have been significant advances in the execution of the binary neutron star inspiral 
MCMC routine since the publication of our previous method jHOJ ; presented here is a 
description of our routine, implemented in C, that employs advanced MCMC techniques 
in order to increase efficiency and performance. 

The outline of the paper is as follows. Section [21 describes our analysis strategy and 
contains a description of Bayesian parameter estimation methods, and the techniques 
we employ in the binary neutron star inspiral MCMC application. In Section 01 we 
present an example application of our method using fictitious data; a binary neutron 
star inspiral event is embedded within noise (that matches the LIGO design sensitivity). 
We conclude the paper with a summary of our results, and a description of our intentions 
to expand this work to more complex situations. 
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2. Analysis strategy 

2.1. Bayesian modeling and MCMC posterior simulation 

Presented here is a short introduction to Bayesian inference, which is quite a different 
way to approach problems like parameter estimation compared to more commonly used 
"frequentist" methods, such as maximum likelihood estimation. The model setup in- 
cludes the specification of a prior distribution for all involved parameters as well as 
the sampling distribution for observed quantities. Inference is done via the derived 
posterior distribution given the observed data; this allows one to make probabilistic 
statements about parameters of interest by treating them as random variables. Due to 
its setup the concept is sometimes dubbed full probability modeling. While the theory 
behind Bayesian modeling is straightforward, its practical application for inference often 
requires the use of numerical methods. 

The sampling distribution is specified in terms of the corresponding probability 
density function (PDF) p(y\9), defining how the observations y come about given a 
(fixed) parameter setting 9. Viewed as a function of 8, it also defines the likelihood 
function. Many frequentist techniques, such as maximum likelihood estimation, are also 
based on the likelihood function, which may be thought of it as 'linking observations to 
parameters'. 

In the frequentist approach 9 is viewed as fixed but unknown. In the Bayesian 
paradigm, on the other hand, 9 is treated as a random variable with a probability 
distribution that reflects the researcher's uncertainty about the parameters. In addition 
to the likelihood a prior PDF p(9) is used, reflecting the pre-experimental knowledge 
and uncertainty about the parameters. The availability of both prior PDF p(9) and 
likelihood p{y\9) then allows to derive the posterior PDF of 9 conditional on the observed 
data y 



by applying Bayes' theorem [31]. Note that the denominator p(y) = Jp(y,9)d9 = 
f p{y\9)p{9)d9 is a normalising constant that does not depend on 9, so the posterior 
PDF p(9\y) is essentially proportional to the product of prior PDF and likelihood 
function. 

The posterior PDF expresses the knowledge about the parameter (s) given model 
(prior PDF and likelihood function) and data. Inference then aims at deriving 
probabilistic statements about parameters (or related quantities) from the posterior; 
one might for example be interested in posterior expectations, quantiles or marginal 
PDFs of quantities of interest, or probabilities of certain events. 

Note that the impact of the prior specification is limited, since asymptotic theory 
shows that the importance of the prior disribution diminishes as the sample size 
increases, so the posterior is eventually dominated by the likelihood fiT\ 132] . When 
MCMC methods are used for posterior analysis as in this paper, the results (MCMC 
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samples) yielded using a certain prior may easily be recycled to generate posterior 
samples assuming different priors by using importance sampling |31j . 

Frequentist and Bayesian analyses are fundamentally different in their approach 
to given problems; Finn jSHj puts the difference as "guessing nature's state" in case of 
frequentist analysis as opposed to "learning from observation" for Bayesian analysis. 
For more detailed discussions see e.g. j2*Tl 1331 E£f] . 

For posterior computations we use a simulation-based approach as the posterior 
distribution is high- dimensional and the calculation of marginal summaries like marginal 
posterior means would require high-dimensional integration. The simulation-based 
approach is to generate a sample from the posterior and then approximate the desired 
integrals by sample averages. For example, the expectation of a random variable X 
with PDF px(x), E[X] = / xpx(x) dx can be approximated through a sample average 
x = j; YliLi x i where x\, . . . , xn are random draws from px(x). Analogously, marginal 
PDFs, quantiles etc. can be estimated from samples from the distribution of interest. 

Markov chain Monte Carlo (MCMC) methods in general may be used to generate 
random sequences of numbers that have a specified stationary distribution and in which 
each random step in the sequence only depends on the previous stage of the chain. 
The Metropolis-algorithm in particular can be applied if the target PDF is known 
only up to a normalising constant. Note that this is the case in equation as the 
evaluation of the numerator p(y) would involve high-dimensional integration. While a 
Metropolis sampler is guaranteed to function as the number of samples goes towards 
infinity, its proper behaviour and efficiency within finite time heavily depends on the 
sensible specifications of its starting point and proposal distribution [3T] . 

2.2. Likelihood, signal templates & priors 

The observed data are a (real-valued) time series of length N and sampling rate ^- 
indicating the phase shift between the two interferometer arms over time at evenly 
spaced time points. In order to analyse a data set, it is Fourier-transformed, from 

{z(t) G E : t = 0, A t , 2A t , ...,{N- 1)AJ (2) 

into 

{2(/)6C: / = 0,A / ,2A / ,...,(f -l)A f }, (3) 
where Af = j^- and 

N-l 

z(f) =A t Y, <3A t ) exp(-27ry7). (4) 

3=0 

The likelihood function depends on the data and a signal template, describing the theo- 
retically derived detector response s$ for a given parameter set *&. The likelihood of 
some parameter set i? then is proportional to the sum of the squared and normalised 
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differences between Fourier transforms of observed signal (z) and signal template (s#) 
over the discrete set of Fourier frequencies {(j x Aj) : ji < j ' < ju} |33] 155]: 

data template 

t ion / 2 4^ XA;)-'^' X A,)| 2 \ 

pW = yx exp^-— 4L j (5) 

noise PSD 

where j^x Af and ju x Af are the lower and upper bounds of the examined frequency 
range, 5 n (-) is the (estimated) noise power spectral density (PSD), NA t is the length 
of the analysed data segment, and K is a normalising constant. 

The signal templates we use are 2.0-post-Newtonian stationary phase approxima- 
tions of the Fourier-transformed inspiral waveform [37] . The model has five parameters: 
the individual masses of the two involved companions mi and 777,2, the coalescence time t c , 
the coalescence phase O an d the effective distance oIe- The effective distance is not the 
actual distance to the source, but also reflects the effect on the amplitude of the gravity 
wave from other parameters and is in general greater than the actual distance. The 
other parameters affecting the gravity wave amplitude are the binary inspiral system's 
sky position and orbital plane inclination angle, the interferometer orientation, and the 
polarization of the wave. We also ignore any effects due to the spins of the compact 
objects. The template is then defined in terms of total mass m t = mi + m2 and mass 
ratio n = =^ as: 

Mf) = o^*» rl ex P(- i (^(/) + 0o + 2nft c )) (6) 

«E 2V6 7T3 C2 v v « ' J 

phase evolution 

W):=X>Ci(/), ( 7a ) 
i=i 

01 = m"' 1 ' (n) 



where 



a-2 



1 /3715 rr \ i 

3M V {^r + 55 V q ■ (7c) 

1 2 

a 3 = 487rg"3, (7d) 

3 128/7 H ' y J 

3 /15 293 365 27145 3085 



a 4 = 1 ?7 H 77 2 o 3 (7 e ) 

128?] V 508 032 504 ' 72 ' ; y ' 



Ci(/) = /"*, C 2 (/) = f-\ Cs(/) = &(/) = /"*, and q = nGm t c' 3 |33 BZ| ■ 

Priors are specified with respect to some preliminary considerations. The coa- 
lescence time t c is assumed to be known in advance with a certain accuracy through 
preprocessing of the data ^01 CHI EH] ; we set its prior to be uniform across ±5ms around 
the true value (which of course is known for our simulated data). The prior for the 
coalescence phase (fio is uniform across its domain [0,27r]. 
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The prior for masses m 1; m 2 and effective distance oIe is set in order to reflect the 
distribution of parameters of a neutron star inspiral given that it has been detected in 
the first place. Initially, the prior for two companions' individual masses (mi and 777-2) 
is uniform between 0.6 and 3.0 M Q (solar masses: M Q rs 2 x 10 30 kg), which effectively 
covers the range of values expected for binary neutron star systems. The prior for oIe 
is derived from the assumption that inspirals happen uniformly across space, so that 
P{d-E < x) oc x 3 . So far, this leads to an improper PDF (that has an infinite integral). 
For actual data, in order to be investigated at all, an inspiral event needs to emit a signal 
that is strong enough to produce a trigger in an inspiral search routine [TEJ [THl d] , and 
thus cannot originate from arbitrarily great distances. We incorporated this constraint 
by downweighing the low- mass/great-distance inspirals according to their lower detection 
probability. The probability of detection is assumed to depend simply on the signal's 
amplitude, which is affected by the inspiral's masses and distance and is proportional 
to: 

a(mi, m 2 , d E ) = -(log(mi) + log(m 2 )) - - log(mi + ra 2 ) - log(d E ) (8) 

2 6 

5 




(mi + 7772) 6 dE/ \ <^e ) 

(so a denotes the logarithmic amplitude; see also equation ©). We could have set 
a threshold amplitude below which neutron star inspirals would be assumed to be 
undetectable, but favoured a 'smoother' transition that does not explicitly assign zero 
probability to parts of the parameter space. We model the dependence between a given 
(logarithmic) amplitude x and detection probability using a (sigmoidal) logistic function 
of the form: 

d a,b(x) = — 1 , x _ a . (9) 

1 + exp(^) 

where a and b are set so that d a ,b(xu) = 1 — p and d a ^(xL) = p for some < p < 0.5 
(e.g.: p := 0.1) and upper and lower reference points xy and xl- So x v denotes the 
amplitude at which the detection probability reaches (1 —p), and x L is the amplitude 
where the probability falls below p. In order to fit d through these points, its parameters 
are set to: 

XL + Xu Xu - x l 

-==-5- and b^J^Ty (10) 

So, given the above prior occurrence and detection probabilities, the resulting (proper) 
joint prior PDF for individual masses m\ and m 2 and distance &e is: 

p(m ll m 2 ,d E ) oc I[ . 6 ,3.o](^i) x I[o.6,3.o] (^2) x d% x d a>b (a(m 1 ,m 2 ,d E ))(ll) 

where 1.4 (x) := I ^ *^ is the indicator function for any set A. For the examples 

[ otherwise 

given in this paper we have chosen values for the sigmoidal logistic function so that a 
2-2M inspiral would be detectable out to 90 and 95 Mpc with probabilities of 90% and 
10% respectively. These values effectively produce a prior PDF for the effective distance 
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that smoothly covers the values of interest, and then falls to zero as the distance gets 
too large. 

Fig. ^ shows estimates of the marginal prior PDFs of masses and effective distance, 
in terms of individual masses (m 1; m 2 ) as well as in the alternative parametrisation of 
chirp mass (m c = ^^\ , 2 ) and mass ratio (77 = t^^^w)- Greater individual masses 



individual masses (m,, m 2 ) effective distance (d E ) 




0.5 1.0 1.5 2.0 2.5 3.0 50 100 150 



chirp mass (m ) mass ratio (n) 




0.5 1.0 1.5 2.0 2.5 0.14 0.16 0.18 0.20 0.22 0.24 



Figure 1. Marginal prior PDFs for masses and effective distance. The upper left 
plot shows the PDF of individual masses (mi, TO2), the lower plots show PDFs in the 
alternative mass parametrisation of chirp mass (m c ) and mass ratio (77). 

are more likely to be observed since low-mass inspirals need to be close to be detectable, 
while high-mass inspirals may as well originate from greater distances: a feature that 
is also known as the Malmquist effect [SB]- Note that the marginal prior PDF for oIe 
extends well beyond 95 Mpc, because the thresholds xl and %u were specified for 2-2M 
binary neutron star inspirals; inspirals involving greater masses still generate signals of 
observable amplitudes when happening further out. 

2.3. Technical details 

The MCMC sampler is implemented in C. Data is imported from the Frame File format 
using the Frame Library jSHj- Before transformation, the data are low-pass-filtered and 
downsampled (from 16 384 to 4096 Hz) |40j : required Filter coefficients are determined 
by the 'Parks-McClellan'-algorithm (see e.g. jS]), again using existing software The 
Discrete Fourier Transform of the data (for both likelihood computation and spectral 
density estimation) is carried out using the free FFTW library 03] ■ The noise spectrum 
is estimated from a data section that is disjoint from the actually investigated data |44j . 
(Pseudo-) random numbers are generated using the Randlib library [4*3] . 

The frequency range over which the likelihood is computed (see equation (jSJ)) was 
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set to 40-1800 Hz. In the original mass parametrisation of individual masses m\ and m-i 
the two masses showed high posterior correlation, making sampling from the posterior 
difficult. In order to improve sampling we expressed the likelihood in terms of chirp 
mass (m c ) and mass ratio (77) instead. In addition we reparameterised the effective 
distance from to \og{d,E), thus implicitly yielding an unbounded parameter space 
and proposal step widths that are proportional to the distance (dg) itself. 

Approximate posterior samples as starting values for the MCMC chains are 
generated using importance resampling, by first generating a large sample from a 
distribution covering the whole prior, and then drawing the actual sample out of these 
with correspondingly assigned weights depending on the posterior PDF [HUES]- The 
Metropolis sampler's proposal distribution was chosen to be a Multivariate Student-t 
distribution. The t-distribution has a similar shape to the Normal distribution, but 
posesses an additional degrees-of-freedom parameter u, and for v — > 00 approaches a 
Normal distribution. It has 'heavier tails' than a Normal distribution, which means that 
extreme values are more likely under a t-distribution than they were under a Normal, 
making it a more robust choice as a proposal distribution [31] • The degrees of freedom 
were set to v := 3, the lowest possible integer value for which the distribution's variance 
is finite. Starting off from an initial setting, the covariance parameter of the proposal 
distribution is recursively adapted to the sample covariance of generated samples during 
an initial burn-in phase [HI EH] . The scale of the proposal covariance is set to ^ of the 
sample covariance estimate, yielding a reasonable acceptance rate for the sampler. 

Convergence of the Markov chains is monitored by using several chains that are 
run simultaneously from different starting points, so one may e.g. verify whether these 
eventually end up in the same mode jlH]. In order to minimise the effect of correlated 
draws from the MCMC output, only every 50th MCMC draw is written to a file which 
is then imported into R, a statistical software, for eventual analysis j^Dl- The C code 
generates some 2500 MCMC samples per minute on a 3.2 GHz Pentium 4 desktop PC. 
The univariate marginal PDFs that are shown in this paper are kernel density estimates 
(for more details about these see e.g. [HI])- 

3. Example application with simulated data 

We illustrate the results of a run of the MCMC sampler on simulated data for which 
the true parameter values are known. The signal analysed had an effective distance of 
25 Mpc, and was embedded in Gaussian and stationary noise that had its noise power 
spectral density match that of LIGO's target sensitivity [22]- The embedded signal 
had a signal to noise ratio of 10. Six parallel chains were run; the starting points of the 
chains were generated by importance resampling of 100 000 draws, a number that proved 
to yield enough eventual draws that were sufficiently close to the main posterior mode 
to ensure reliable and fast convergence of the Metropolis algorithm. The first 30 000 
iterations of each chain were considered the burn-in-phase, during which the iterations 
15 000-30 000 were used to tune the proposal covariance. The code was then run for 
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2 million iterations in total, which after thinning out of the samples and discarding the 
burn-in yielded a sample of size 236 400 from the posterior. The 'multivariate potential 
scale reduction factor' BP was close to 1 (B p = 1.0034), indicating convergence of all 
chains [55] . 

Figures EHSl illustrate estimates of marginal posterior PDFs, that is, univariate 
PDFs with respect to one specific parameter, with all other parameters integrated out. 
Firstly, Fig. |21 shows posterior PDFs for the five individual parameters. Most of them 
exhibit a mode near the true parameter value (indicated by dashed lines). One can see 



chirp mass (m c ) mass ratio (11) 




1.212 1.214 1.216 1.218 1.220 0.18 0.20 0.22 0.24 



coalescence time (t c ) coalescence phase (4> ) effective distance (d E ) 




84.614 84.616 84.618 84.620 84.622 JI 2lt 20 25 30 35 40 45 



Figure 2. Marginal posterior PDFs of the five parameters. Dashed lines indicate true 
parameters. 

that the relative precision of parameter estimation varies significantly between different 
parameters. For example, the posterior of the chirp mass covers a range of about 
0.006 M Q , while the prior range initially was some 2.1 M Q (cp. Fig. Q). The coalescence 
phase's posterior, on the other hand, still covers the complete prior domain. 

Fig. 01 allows for some insight into joint PDFs of some of the parameters. The 
joint PDF of chirp mass and mass ratio (j5k) shows a positive correlation between the 
two parameters. Fig. shows interaction between two parameters (</>o and i]), and in 
particular demonstrates that although the marginal PDF of <po alone is almost uniform 
(see Fig. |2J), this does not imply that its effect on the posterior was negligible. 

The MCMC sampler internally works with chirp mass (m c ) and mass ratio (r/) 
instead of individual masses (mi, m 2 ). A posterior sample of the individual masses still 
can easily be obtained by back-transforming each pair of (m c ,r]) samples. Fig. 0] shows 
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a) chirp mass (m c ) vs. mass ratio (r|) b) coalescence phase (<M vs. mass ratio (r|) 




1.210 1.212 1.214 1.216 1.218 it 2it 

Figure 3. Bivariate marginal posterior PDFs for two pairs of parameters. Dashed 
lines indicate the true values, the true coalescence phase is (f>o — (Histograms, the 
greyscale plots show relative PDFs normalised to the mode). 



these two marginal PDFs combined into one plot. 




1.0 1.5 2.0 2.5 



Figure 4. Combined plot of marginal posterior PDFs of the two companions' indi- 
vidual masses (assuming mi < 771.2). Dashed lines indicate true parameters. 

Analogously, other functions of the parameters can be derived and distributional 
features investigated; if e.g. one was interested in whether the masses differ 'significant- 
ly' or are 'almost equal', we can estimate: P(m 2 > 3m!) = 0.11% or P(m 2 < 1.5mi) = 
4.84%. Fig. El shows the posterior PDF of the logarithmic amplitude a(mi,77i2, djs) 
(JHJ). Comparing it to the prior PDF you can see that, since it is significantly above 
the reference points xy and x^, the particular specification of the lower bound of the 
parameter space does not affect our conclusions. Table [T] shows summary statistics of 
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l = a(2, 2, 95) 

Xy = a(2, 2, 90) / 


a(1 , 2, 25) 













-4.0 -3.8 -3.6 -3.4 -3.2 -3.0 -2.8 



Figure 5. Posterior PDF of the signal's logarithmic amplitude a(mi,m2, de); the 
dashed line indicates the true value. The prior PDF (dotted line) and xjj and xl are 
shown as well. 



the posterior distributions of the inspiral's parameters. 

Table 1. Posterior estimates: Means, medians and 95% central posterior intervals for 
several parameters. 



parameter 


mean 


median 


95% c.p.i. 


true 


unit 


chirp mass 




1.2161 


1.2159 


[1.2145,1.2186] 


1.2167 


M 


mass ratio 


(v) 


0.2174 


0.2162 


[0.1987,0.2457] 


0.2222 




coalescence time 


(tc) 


84.6174 


84.6174 


[84.6160,84.6189] 


84.6167 


s 


coalescence phase 






- not meaningful - 


0.0 


radian 


effective distance 


(d E ) 


26.28 


25.99 


[21.55,32.68] 


25.00 


Mpc 


mass 1 


(mi) 


0.980 


0.964 


[0.876,1.229] 


1.0 


M 


mass 2 


("i 2 ) 


2.062 


2.085 


[1.600,2.327] 


2.0 


M 



4. Discussion 

We have developed an efficient method for examining gravitational wave interferometer 
data for binary neutron star inspiral signals. Our MCMC code investigates the five 
parameters of an inspiral signal that is described by its two masses, effective distance, 
time, and phase at coalescence. This code can be applied at the end of a binary 
neutron star detection pipeline [To! HHj IT^j . thereby extracting detailed information 
about parameters and related quantities. 

We are extending our MCMC research on binary neutron star inspiral signals. 
The power of MCMC methods is their ability to address complex signals with large 
parameter numbers. We see MCMC methods being of great use for even more complex 
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gravitational wave problems; an example would be inspiral signals from black hole-black 
hole, or black hole-neutron star systems, and binaries with significant amounts of spin 
for each mass. In order to address these signals higher order post-Newtonian terms are 
necessary |54j . A current and promising research effort of our group is in incorporating 
these high order PN terms, and to create an efficient MCMC search routine that will 
find all of the signal parameters. Another on-going research project is in applying 
MCMC methods to the multiple detector binary neutron star inspiral problem, where 
sky position parameters can also be estimated. We expect that MCMC methods will 
prove useful with more complex binary inspiral scenarios. 
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